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ABSTRACT 

Yttria-stabilized zirconia (YSZ) is of interest to the aerospace community, notably for its application 
as a thermal barrier coating for turbine engine components. In such an application, diffusion of both 
oxygen ions and cations is of concern. Oxygen diffusion can lead to deterioration of a coated part, 
and often necessitates an environmental barrier coating. Cation diffusion in YSZ is much slower 
than oxygen diffusion. However, such diffusion is a mechanism by which creep takes place, 
potentially affecting the mechanical integrity and phase stability of the coating. In other applications, 
the high oxygen diffusivity of YSZ is useful, and makes the material of interest for use as a solid- 
state electrolyte in fuel cells. 

The kinetic Monte Carlo (kMC) method offers a number of advantages compared with the more 
widely known molecular dynamics simulation method. In particular, kMC is much more efficient for 
the study of processes, such as diffusion, that involve infrequent events. 

We describe the results of kinetic Monte Carlo computer simulations of oxygen and cation diffusion 
in YSZ. Using diffusive energy barriers from ab initio calculations and from the literature, we 
present results on the temperature dependence of oxygen and cation diffusivity, and on the 
dependence of the diffusivities on yttria concentration and oxygen sublattice vacancy concentration. 
We also present results of the effect on diffusivity of oxygen vacancies in the vicinity of the barrier 
cations that determine the oxygen diffusion energy barriers. 
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Motivation 


Yttria-stabilized zirconia (YSZ) is of interest to the aerospace community, notably for 
its application as a thermal barrier coating for turbine engine components. In such an 
application, diffusion of both oxygen ions and cations is of concern. Oxygen diffusion 
can lead to deterioration of a coated part, and often necessitates an additional anti- 
oxidation coating. On the other hand, the high oxygen diffusivity of YSZ makes the 
material of interest for use as a solid-state electrolyte in fuel cells. Cation diffusion in 
YSZ is much slower than oxygen diffusion. However, such diffusion is a mechanism 
by which creep takes place, potentially affecting mechanical integrity and phase 
stability. 


Why Kinetic Monte Carlo? 


Currently, molecular dynamics (MD) is the method of choice when examining the 
dynamical behavior of many-atom systems. MD produces detailed trajectories of 
atoms, but is typically limited to time scales on the order of nanoseconds.When the 
phenomenon under study involves so-called "infrequent events," MD may be very 
inefficient. A diffusive hop (the mechanism behind both oxygen ion and cation 
diffusion) is such an event.kMC is a simulation method designed to deal with 
infrequent-event systems. It averages out the system's behavior during 
"uninteresting" times and concentrates on the events of interest, and is typically 
much faster than MD for such systems.lt requires a knowledge of all the "interesting 
events accessible to the system.lt does not provide detailed atomic trajectories. 


We describe the results of kinetic Monte Carlo computer simulations of oxygen and 
cation diffusion in YSZ. Using diffusive energy barriers from ab initio calculations 
and from the literature, we present results on the temperature dependence of 
oxygen and cation diffusivity, and on the dependence of the oxygen diffusivity on Y 
concentration. We also present preliminary results of the effect on diffusivity of 
oxygen vacancy clustering near Y ions. 



kMC, as applied to diffusion 

The sole events of interest are diffusive hops among vacancy sites on the oxygen or cation 
sublattices. 

Because the vibrational frequency is much larger than the average frequency of hops, 
the system loses its memory of the previous event (here, which of a vacancy’s neighbors 
were involved in the previous hop). Therefore, the events can be considered as 
independent. 

In such cases, the probability per unit time that a vacancy will undergo a diffusive 
hop is constant, and the probability distribution for the first hop is given by 

p{t) = v to texp(-v tot t) 

where v tot is the aggregate rate constant for all hops accessible to the vacancy. 

Because all events are independent, the overall hopping rate is the sum of the 
individual rates for each accessible hop. 

In diffusive hopping, the rates depend on the energy barriers: 

vab = v°exp(-E AB /k B T) 

in which u ab and Eab are the hopping rate and migration barrier energy for a 
hop between oxygen or cation sites A and B respectively, and is the frequency 
factor. 

The hopping probability is computed from the hopping rat £*ab - 
where r is the sum of all hopping rates. 


A computational cell is constructed with the appropriate numbers of Y 
cations and O vacancies distributed randomly, consistent with charge 
neutrality. 

All events accessible to the system are identified and the corresponding event 
probabilities are computed and collected in an event catalog. 

Iterative kMC procedure: 

An event is chosen at random and executed. 

The simulation clock is advanced by a time stdpfc = — 

where R is a random number greater than zero and less than or equal to unity. 

The event probabilities and rates are modified to reflect the new state of the 
system, and are included in the modified catalog. 

When the simulation has run long enough to accumulate statistically useful information, 
the mean square displacement, averaged over all vacancies, is computed. The vacancy 
diffusivity is obtained from the mean square hopping distance using the Einstein v 
relation 

< R 2 >— 6D v t 

and the ionic diffusivity Di is obtained by balancing the number of vacancy and ion hops: 

n . — n 

— i _r* 

-L L-H; 

where C lv is the concentration of vacancies on the sublattice of the hopping ion. 



Oxygen diffusivity in YSZ 

Oxygen diffusivity of YSZ over the temperature range of interest is typically in the 
range of 1 0" 6 - 1 0' 7 cm 2 /sec. 

Diffusivity increases with increasing temperature. 

Diffusivity increases with increasing Y concentration at low concentrations, and reaches 
maximum at 8-12 mol percent Y 2 O 3 . 

Results from experiment 

Oishi and Ando find oxygen diffusivities of 0. 5-3.0 x 10' 6 cm 2 /sec over a 
temperature range of 1373K-2073K. 

Strickler and Carlson find diffusivities in the range of 0.4-4. 6 x 10' 6 cm 2 /sec, and 
observe a maximum in the diffusivity as a function of concentration. They find that 
the maximum occurs at slightly higher concentrations at higher temperatures. 

Results from simulation 

Khan et al., produce diffusivities of 1 .9 x 10' 7 to 3.0 x 10' 6 cm 2 /sec over the 
temperature range 873K-2073K. 

Krishnamurthy et al., produce diffusivities of 0.1 -7.0 x 10' 6 cm 2 /sec over a 
temperature range of 1000-2200K. 
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Diffusion path is to the nearest neighbor on oxygen sublattice (100 direction). Each 
site on the oxygen sublattice has four cation neighbors that form a tetragon, and 
hopping occurs across tetragonal edges, with the hopping oxygen ion passing 
between two barrier cations. 

Barrier energies are assumed to depend only on the two barrier cations. The possible 
barrier pairs are Zr-Zr, Zr-Y and Y-Y. 



Oxygen ion hopping path. Oxygen ions are green, cations grey. 

MD simulations have yielded activation energies of 0.37 eV (Khan) and 0.2-0. 8 eV (Li 
and Hafskjold). A tracer diffusion study finds a value of 0.44 eV, while results from bulk 
conductivity and ac impedance spectroscopy give values of 0.79-1 .12 eV. 


X. Li and B. Hafskjold, J. Phys. Condens Matter 7, 1255 (1995) 




Energy barriers computed from ab initio Car-Parinello MD by Krishnamurthy et al. 
are shown below, along with values used in this work, obtained using the Abinit 
plane wave pseudopotential code with pseudopotentials from the Fritz-Haber 
Institute. Note that the ordering of the barrier energies is Zr-Zr, Zr-Y and Y-Y. 


Ab Initio Oxygen Diffusion Energy Barriers (Energies in eV) 


Authors 

Method 

Zr-Zr 

Zr-Y 

Y-Y 

Krishnamurthy, 
Yoon, Srolovitz 
and Car 

Car-Parinello 

MD 

0.58 

1.29 

1.86 

Krishnamurthy, 

Srolovitz, 

Kudin and Car 

Car-Parinello 

MD 

0.473 

1.314 

2.017 

This work 

DFT 

0.706 

1.214 

1.941 



















kMC diffusivities as a function of Yttria concentration are shown below. The absolute values 
of the diffusivity are in general agreement with experiment and with other simulations. 
Notably, the trend is consistent with experiment, and the maxima are in approximately the 
right places, and the location of the maximum increases with increasing Y concentration, as 
observed experimentally. 


This trend has been postulated to be 
the result of a competition between two 
effects. 

(1) As Y concentration increases, the 
number of oxygen vacancies increases 
so as to maintain electrical neutrality. 
Increasing the available vacancies 
increases the diffusivity.(2) As Y 
concentration increases, the number of 
relatively high-energy Zr-Y and Y-Y 
barriers increases, inhibiting diffusion. 


Oxygen Diffusivity Versus Y Ion Concentration 
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Yttria Molar Concentration 


Temperature dependence of the oxygen diffusivity for a number of Yttria 
concentrations is shown below. Activation energies range from 0.25eV to 0.46eV, in 
good agreement with experiment. 

The activation energies increase with increasing Yttria concentration. This is 
consistent with the larger numbers of high energy Zr-Y and Y-Y barriers present at 
higher Yttria concentrations. 

Temperature Dependence of Oxygen Diffusivity 
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There is some indication that oxygen vacancies may preferentially occupy either 
yttrium or zirconium nearest neighbor sites. If so, this may affect the barrier 
energies. We have computed the change in barrier energy for two cases. 

(1) A yttrium ion is a nearest neighbor to a hopping oxygen ion. 

(2) An oxygen vacancy is a nearest neighbor to a yttrium barrier ion. 


Barrier energies with preferential O vacancy-Y ion clustering (energies in eV) 



Zr-Zr 

Zr-Y 

Y-Y 

Original 

0.706 

1.214 

1.86 

Y adjacent to 
hopping 0 

0.691 

1.140 

1.941 

0 vacancy 
adjacent to Y 
barrier ion 


1 .0006 

1.64 

















There appears to be only a small effect 
on diffusivity at large Y concentrations, 
but additional configurations are currently 
being investigated. 


Diffusivity, cm /sec 


Oxygen Diffusivity versus Y concentration, 
Effects of Y-0 Vacancy Clustering 
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Cation diffusion 

Cation diffusion is treated in a similar way, except that the diffusion path and 
energetics are different. 

Initial ab initio calculations, along with MD results from the literature, suggest that nearest 
neighbor hops in the (110) direction are energetically favorable compared with longer 
(100) hops. 

The barrier of oxygen ions consists of a single pair. 



Cation hopping path. Oxygen ions are red, cations grey. 




Ab initio calculations to date yield barrier energies that are somewhat too large. A more 
detailed ab initio reaction path calculation may resolve the issue. 

Experimental values of the hopping enthalpies, and the ordering of the Zr and Y 
enthalpies, are not consistent. 

We have chosen barrier 
energies of 3.7 and 3.62 for 
Y and Zr migration, 
respectively. These values 
have been chosen to 
produce diffusivities within 
the range of experimental 
values at 2000K. 

Note that the difference in 
barrier energies for Zr and Y 
hops is small. 
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Author 

Temperature 
Range, C 

Value, eV 

Solmon 

1300-1700 

4.8-4.95 

Gomez-Garcia 

above 1500 

5. 5-6.0 

Chien and Heuer 

1100-1300 

5.3 

Dimos and Kohlstedt 

1400-1600 

5.85 

Kilo 

1125-1460 

4. 4-4. 8 

Mackrodt 


2-7eV (migration 
energies) 

This work 


3.7 (Y), 3.62 (Zr) 




Diffusivity is expected to depend on cation vacancy concentration and on Y 
concentration. 


Cation diffusivity is expected to increase with cation vacancy concentration, but the number 
of such vacancies is small. 


Because the barrier energies for Zr and Y hopping are not very different, the direct 
effect of increasing the Y concentration is expected to be small. 

Increasing the Y concentration increases the oxygen vacancy concentration, which 
increases the probability that a vacancy will exist in the barrier oxygen pair. Ab initio 
calculations suggest that such vacancies will reduce the barrier energy, and increase 
the diffusivity. 

Results from experiment and simulation are not conclusive, but only a weak 
concentration dependence is typically observed. 



The vacancy concentration on the cation sublattice is much smaller than on the oxygen 
sublattice. We consider cation vacancy concentrations of one percent or smaller. 


Diffusivity versus temperature for our 
results and results from the literature 
are shown at right. The slope of the 
line from our kMC calculations is less 
than that of the experimental results 
shown, although it differs from the 
results of Chien et al. by less than ten 
percent. Other kMC simulation results 
(not shown) using larger barrier 
energies show a slope more 
consistent with other experimental 
results the reasonable agreement 
between simulation and experiment of 
the slopes of the ln(D) versus 1/T 
suggests that the kMC simulations 
capture the fundamentals of the 
diffusion process. 


Temperature Dependence of Cation Diffusivity 



Conclusions and Outlook 


kMC is capable of producing accurate simulations of "infrequent-event" phenomena, 
and provides insight into dynamical behavior without detailed atomic trajectories. 

Oxygen diffusivities from kMC simulations using ab initio barrier energies are in 
reasonable agreement with experiment and other theoretical work. 

There are small differences in the barrier energies depending on the cation species 
adjacent to the hopping oxygen ion or the barrier cations. Therefore it may be 
problematic to assume that oxygen vacancies and dopant cations are randomly 
distributed within the model lattice. It may be possible to generate more realistic 
distributions using a Metropolis Monte Carlo species exchange algorithm. 

Ab initio barrier cation hopping energies for the assumed barrier complex are too large. 
Reaction path calculations may provide insight. 

Using barrier energies derived from experiment produces cation diffusivities and 
temperature dependences that are in reasonable agreement with experiment. 

We are currently extending our work to other dopants, and other TBC and EBC 
candidate materials. 

Enlightening discussions with Beth Opila, Dongming Zhu, Nate Jacobson and Jim 
Smialek are gratefully acknowledged. 


Possible additions to conclusion slide 


KMC simulations can produce realistic diffusivities and temperature and 
concentration dependences. 



kMC, briefly 


All possible events of interest accessible to the system are identified, 
and the energetics computed. For each event, a probability, and the 
associated rate, are computed.AII events and probabilities are organized 
into a catalog. Iteratively: An event is chosen at random from the 

catalog and executed. The simulation clock is advanced 
stochastically, consistent with the rates. 

The event catalog is updated to reflect the new accessible events and 
their probabilities. 






